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Dynamical models of the elliptical galaxy NGC 4494 
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ABSTRACT 

We present dynamical models of NGC 4494, which we built using our iterative method 
presented in a previous paper. These models are live ./V-body models consisting of equal 
mass particles, and they are steady state as confirmed by a fully self-consistent evo- 
lution. Our goals were twofold. The first one - namely to test whether our iterative 
method could indeed be used to construct galactic models following given observa- 
tional constraints, both photometric and kinematic - was fully achieved. Our method 
allowed us to go beyond a simple spherical model and to make full sets of rotating, 
axisymmetric models without any limitations to the velocity distribution. Our second 
goal was to understand better the structure of NGC 4494, and more specifically to 
set constraints on its halo mass. For this we tried three families of models: without 
halo, with a light halo and with a heavy halo, respectively. Our models reproduce 
well the photometry and the kinematics, the latter except specific regions where some 
non-equilibrium or non-axisymmetric structure could be present in the galaxy (e.g. 
the kinematically decoupled core). However, the lower order moments of the velocity 
distribution (up to and including the second order) do not allow us to discriminate 
between the three halos. On the other hand, when we extend the comparison to the 
higher order moments of the velocity distribution obtained from the long-slit data, we 
find that our light halo model fits the data better than the no halo, or the heavy halo 
models. They also reproduce the shape of the angular dependence of the PNe velocity 
dispersion in the outermost parts of the galaxy, but not the amplitude of its azimuthal 
variation. This may imply that a yet more general class of models, such as triaxial, 
may be necessary for a yet better fit. 

Key words: methods: numerical - methods: N-body simulations - galaxies: elliptical 
and lenticular - galaxies: individual: NGC 4494 - galaxies: kinematics and dynamics 



1 INTRODUCTION 

Dark matter around ordinary elliptical galaxies is one of 
the hottest topics in dark matter studies today. The main 
goal is to obtain sufficient constraints on the dark matter 
mass from observed stellar kinematics. Traditional long-slit 
absorption line spectroscopy can only very rarely give kine- 
matics outside 2R e , where R e is the effective radius encom- 
passing half the tot al light of the galaxy (see for example, 
ICoccato et al.l l2009h It is, nevertheless, possible to obtain 
line-of-sight velocities at larger radii using planetary nebulae 
(PNe), because their strong emission line at 5007 A [OIII] 
stands out against the faint galaxy background. It is usually 
assumed that PNe trace the kinematics of the underlying 



field star^j. It is thus possible to obtain from the PNe the 
stellar kinematic parameters at the per i phery of the galaxy, 
out to 5 — 7 R e JGoudfrooii et al.ll 1994: R omanowskv et al.l 



20031: iDouglas et al.l 20071: Ide Lorenzi et all 120081 , 120091 ; 
Napolitano et al . 2009; C occato et a,l.ll2009l ). 

Romanowskv et al.l (|2003l ) studied the kinematics in the 



outer part of three ordinary elliptical galaxies, namely NGC 
821, 3379 and 4494 (out to 4 - 6 R e ) and found that their 
velocity dispersion profiles decline nearly Keplerian-like at 
radii outside 2R e . They modeled the observational data by 
means of sp herical Jeans models and b y means of orbit based 
models ( see iRomanowskv et al.l IMPOST ) for details) and they 
showed that their data are consistent only with models with 
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1 iDekel et al.l J2005T) , however, noted that observations of PNe 
can be biased towards the kinematics of younger stars. 
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little or no halo. This result, in go od agreement with what 
was already found for NGC 3379 bylGoudfrooii et al.l || 19941 ) 
and for NGC 4697 bv lMendez et al.l (J200ll ). is vervsurpris- 
ing. Indeed, present theoretical and observational data argue 
that ellipticals are formed from mergings of spirals, which 
are known to have a considerable amount of dark matter 
i|Bosmall2004l . and references therein). So, if the progenitors 
have a considerable amount of dark matter, how can the 
merger product not have it? Furthermore, this result is in 
conflict with the predictions of the standard ACDM cosmol- 
ogy- 

iDekel et al.l (|2005l ) constructed elliptical galaxy models 
from numerical simulations of mergers of a pair of disk galax- 
ies. Their resulting models have a "normal" massive dark 
halo and a velocity distribution with a high radial anisotropy 
in the outer parts. The latter leads to a low observed (from 
most viewing angles) velocity dispersion in the outer parts, 
which in turn leads to a low estimated halo mass, contrary to 
the r eal dark halo mass of the model, which is normal. In this 
way, Dekel et al . (2005) e xplain the resu lts of iMendez et al.l 
(J2001) and of iRomanowskv et al.l 1120031 ) as due to the ve- 
locity distribution in the merger remnant. 

Since thi s velocity anisotropy is so crucial to the data in- 
terpretation, lAthanassoulal (J2005J) examined whether it was 
general, or whether it depended on the specific mergers used. 
She examined the velocity anisotropy in multiple mergers, 
as would occur e.g. in groups. In such cases, a pair merger is 
not examined as an isolated event, but a whole sequence 
of mergers is considered. This model would be more re- 
alistic in groups, but is also in good agreement with the 
standard ACDM cosmology. The result of such mergers are 
com patible with the observed properties of elliptical galax- 
ies (IWeil fc Hernquistlll994l . Il99a lAthanassoula fc Vozikisl 
1999) . Concerning the velocity anisotropy in the outer parts, 
she found that the result is more complex than the single 
pair mergers would predict and that the anisotropy depends 
strongly on the time between two successive mergers. Thus, 
m ore work is necessa ry to establish how general the result 
of IDekel et al.l (120051) is. A similar conclusion was reached 
bv iDouelas et all (J2007) who modeled the data of NGC 
3379 and argued that there are considerable discrepancies 
between the observations and dark-matter-dominated sim- 
ulations and re-iterate the question of whether NGC 3379 
has the kind of dark halo that the current ACDM paradigm 
requires. 

Ide Lorenzi et al.l IJ2008I . 120091 ) constructed dynami- 



cal models of NGC 4697 and 3379 



using 



the 



X 



made-to-measure part icle method (|Sver fc Tremaineill996l : 
Ide Lorenzi et al.ll2007h implemented in the NMAGIC code. 
Their main result is that the observational data are consis- 
tent with a fairly wide range of halo mass profiles, although 
it was poss ible to pl ace some limits on the halo mass. For 
NGC 4697, Ide Lorenzi et all (2008) found that models with 
a low density halo with v c (5R e ) < 200 kms -1 are not con- 
sistent with the data, where v c (5R e ) is the total circular 
velocity at 5R e . This, however, is a rather weak limit be- 
cause eyen_t4ie2r_jnodel Dwith v c (5R e ) w 210 (see figure 
15 in Ide Lorenzi et al.l IJ2008I )) has a very light halo which 
contributes only abo ut 35% of the m a ss wi thin 5R e . For 
the galaxy NGC 3379 Ide Lorenzi et all (J2009I ) found that a 
model without a halo, as well as a model with a heavy halo 
with v c {7R e ) > 250 fans -1 would be excluded by the ob- 



servational data, but only at a la confident level. So these 
constraints are not very strong. 

The main problem in determining the halo mass pro- 
file from the observed kinematics in elliptical galaxies is 
the well known mass-anisotropy degeneracy. The low veloc- 
ity dispersion on the periphery of some elliptical galaxies 
can be explained either by nearly isotropic models with a 
light halo or by radially anisotropic models with a heavy 
halo. It is generally accepted that the mass-anisotropy de- 
generacy can be broken by means of high order moments 
of the line-of-sight velocity distribu tion (LOSVD. iGerhardl 
1 19931 ; Ivan der Marel fc FranxHl993r ). As shown in these pa- 
pers, isotropic models have a Gaussian LOSVD and radially 
anisotropic models have centrally peaked LOSVD as well as 
long tails. So these models can be distinguished by means of 
high-order moments of LOSVD. Howe ver, th e hig hly radial 
anisotropic models presented in IDekel et al.l (J2005J) have a 
LOSVD with relatively weak d eviations from Ga ussian (see 
supplementary information in IDekel et al.l 12005). Unfortu- 
nately, this means that, at least in some cases, breaking 
the mass-anisotropy degeneracy can be very difficult, if not 
practically impossible. 

Let us also point out an obvious, but sometimes ig- 
nored, problem. From a mathematical point of view, it is 
possible to prove the existence of a given type of model by 
simply constructing dynamical models, but it is not pos- 
sible to prove its non-existence. Let us, for example, con- 
struct dynamical models of a real galaxy by means of the 
NMAGIC method or by means of our iterative method (see 
below). More specifically, let us construct an axisymmetric 
model with some dark halo. If this satisfies all observational 
data, then we have proven the existence of a model with 
such a halo agreeing with the observational data. But if we 
fail to construct such a model, then formally we have not 
proven anything, since we can not exclude that our failure 
is due to the method itself, or to the fact that we have not 
searched sufficiently. We would have needed to prove that, if 
an equilibrium model with given parameters did exist, then 
our method would construct it. This is not straightforwardly 
proven either for our iterative method, or for the NMAGIC 
method, or for any other orbit based method. Moreover, 
even if we did prove it, we would have only proven that an 
axisymmetric model with this particular halo is excluded by 
the observational data. We would not have proven anything 
about triaxial models, or about models with a somewhat 
different halo mass profile. C onsequently, the conclusion of 
Ide Lorenzi et al.l (|2008l . |2009j) that the observational data 
of NGC 4697 and 3379 are consistent with a wide range of 
halo masses can be firmly believed. But, on the contrary, if 
one finds that a set of models with a heavy dark halo do 
not agree with the observational data, then one has at the 
best only an argument that the galaxies in question have a 
light dark halo, or no halo. Unfortunately, this is no proof, 
since another, more general, type of heavy halo could per- 
haps have fitted the observations. It is thus very useful to 
try different approaches, to see whether this disagreement 
persists, or not. If more than one method lead to the same 
conclusion, then the argument is considerably strengthened. 

Initially the p resen t work was inspired by an article of 
iNapolitano et al.l l|2009l . hereafter N09), where the authors 
presented a large amount of new observational data of the 
ordinary elliptical galaxy NGC 4494, resulting in positions 
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and velocities of 255 planetary nebulae out to seven effective 
radii. They also presented new wide-field surface photometry 
from MMT/Megacam, and long-slit stellar kinematics from 
VLT/FORS2. Using these data they put constraints on the 
distribution of dark matter in this galaxy. They constructed 
spherical dynamical models of the system using two different 
methods, but both of them are based on Jeans equations. 
They argue that some dark matter is required by the data 
and the ir best fit model has a relatively low h a lo ma ss. 

In iRodionov. Athanassoula fc Sotnikoval (120091 . here- 
after RAS09) we presented an iterative method for con- 
structing equilibrium iV-body models with given prop- 
erties. This method has been already widely used to 
construct initial conditi o ns for TV-body simulations (e.| 
IRodionov fc Orlovl 120081 ; iMachado fc Athanassoulal I201C 
and can also be directly used for constructing dynamical 
models of real galaxies from observational data. In contrast 
with NMAGIC models, our models consist of particles with 
equal masses. They are steady-state and can be directly used 
in A^-body simulations. For example, we can directly check 
that the constructed model is indeed in equilibrium. Also, 
we can easily calculate for this model any parameter which 
can be directly obtained from the positions and velocities of 
the particles. 

In this article we will apply our iterative method to 
the construction of dynamical models of NGC 4494 from 
the observational data presented in N09. We have two pur- 
poses. First, our intention is to demonstrate that our itera- 
tive method can indeed be successfully used to construct dy- 
namical models of real galaxies. The second aim is to return 
to the interesting question of the halo mass of NGC 4494, 
using a different method from that of N09, to see whether 
results can in any way depend on the method used. This 
is particularly important since the discrepancy between the 
no-halo model of N09 and the observational data are not 
large (see upper panels of their figure 12). 

We present the observational data in Sect. [2] and our 
method in Sect. [3] In Sect, fj] we describe our models and 
compare them to observations. We summarise and conclude 
in Sect. [5] 



2 PREPARATION OF OBSERVATIONAL 
DATA 

We use the same observational data as N09. These include 
the surface photometry, the stellar kinematics along the ma- 
jor and minor axes obtained by means of long-slit spec- 
troscopy, and velocities and positions of 255 planetary nebu- 
lae (see N09). When constructing our models, we use a phys- 
ical system of units, i.e kpc and Mq. We adopt a distance 
to NGC 4494 of 15.8 Mpc (see N09), so that 1" is equal to 
0.0766 kpc. Our models can be easily rescaled to any other 
distance. Let us assume we have an equilibrium model con- 
structed for some distance di, and we want to rescale it to 
a distance d-z — Cd\. To keep the surface photometry un- 
changed we need to rescale all space coordinates of particles 
as r2 = Cr\. To keep the model in equilibrium we need to 
change the mass of the model as M2 = CM\ . Particle ve- 
locities need not be changed, so all the observed kinematic 
parameters are unchanged. The new total luminosity of the 



2005h. 



HI), 



galaxy is L2 = C 2 L\, so that the new mass-to-light ratio is 
M 2 _ 1 Mi 

Let us now describe how we prepare the observational 
data for use in our iterative method. 



2.1 Surface photometry 

We use the combined V-band surface photometry of NGC 
4494 presented in table Al of N09. Th ese data are a com - 
bination of HST-based observatio ns of Lauer et al.l I 
ground based CCD observations of lGoudfrooii et al.l I 
and the new observations of N09. 

We, furthermore, make the following simplifications: We 
assume that the ellipticity e is the same at all radii and equal 
to 0.162, the mean value found in N09. We also assume that 
the shape of the isophotes is precisely elliptical, so as not to 
introduce in the analysis unconstrained high order isophote 
shape parameters. The first two rows of table Al in N09 give 
the surface brightness as a function of the intermediate axis 
Rm, measured in arcseconds. This is related to the ellipticity 
and to the major axis R a by R m = _R a \/l — e. We convert 
surface brightness in units of Lq^v/pc 2 (assuming an abso- 
lute magnitude of the Sun in the V-band Mq,v = 4.8) and 
R m in parsecs using the adopted distance of 15.8 Mpc. 

Excluding the innermost region (R m < 5"), this s urface 
brigh tness profile is fitted very well by the Sersic law IjSersid 
1968) 



I(Rm) = /oexp(-( J R m /a s ) 1/n ) 



(1) 



with parameters, 7o = 41764 Lq^v/pc 2 , a s — 0.008809 kpc 
and n — 3.3, as shown by N09. 

Numerically, the surface brightness profile is described 
as follows: Inside the region R m < 0.46 kpc (ss 6") we inter- 
polate the tabular data linearly. In the region 0.46 kpc < 
Rm < 40 kpc we adopt a Sersic profile. In the region 
40 kpc < R m < 50 kpc we trunca te the Sersic pr ofile by 
means of a 5-th order polynomial IjDehnenl l2000bl . eq. 4). 
From the profile of the surface brightness and the adopted 
value of ellipticity we can calculate the two-dimensional dis- 
tribution of surface brightness. The total luminosity of this 
model is L v = 2.36 x 10 10 L©,v, or M v = -21.13. We as- 
sume that the stellar mass-to-light ratio (M/L) is constant, 
in which case the distribution of surface brightness equals 
the surface mass distribution to within an unknown multi- 
plier M/L. 



2.2 Kinematical data 

2.2.1 Symmetries and system of coordinates 

In order to prepare the kinematical data so that they can 
be used by our iterative method, we need first to define a 
system of coordinates which will be used for the model and 
to assume what symmetries the galaxy has. 

In all the following we will assume that the galaxy 
is axisymmetric. Elliptica l galaxies can well be triaxial 
IjBinnev fc Tremaine 2008). Triaxial models, however, have 
extra free parameters that add complexity to the modeling 
and are beyond the scope of this paper. 

Let us consider a Cartesian (X,Y ,Z) coordinate system 
such that the sky plane coincides with the XZ plane. We 
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choose the Z axis so that it coincides with the minor axis of 
the projected image of the galaxy and then the X axis will 
coincide with the major axiqj. We assume that the rotation 
axis is perpendicular to the X axis and therefore in the YZ 
plane. This rotation axis is denned by its angle a with the 
Z axis. If a = then the rotation axis coincides with the Z 
axis, so the galaxy is edge-on. 

We also assume that the galaxy has a reflection sym- 
metry with respect to the plane of symmetry perpendicular 
to the rotation axis and centered on the center of coordi- 
nates (0, 0, 0). In that case, the observed image of the galaxy 
will have a reflection symmetry about both the X and the 
Z axes. The observed line-of-sight velocity distribution will 
have a reflection symmetry about the X axis. Also the ob- 
served line-of-sight velocity distribution will have a reflection 
symmetry about the Z axis, except for the velocity sign. It 
means that the points (x, z) and (—a;, z) will have the same 
line-of-sight velocity dispersion and that their line-of-sight 
mean velocities will be equal but with opposite sign. The 
velocity distribution in the points (x, z) and (x, —z) will be 
fully identical. Consequently, if we know the line-of-sight ve- 
locity distribution in any of the four quadrants, we know it 
for the whole system. Similarly, if we assume such symme- 
try then all the observed kinematical data can be "reduced" 
to the first quadrant of the sky plane (x > 0, z > 0). So 
in a first stage, we reduce all kinematical data to the first 
quadrant. 

In our iterative method the input kinematical data 
should be given as mean velocities and velocity dispersions 
in a set of two-dimensional areas on the sky plane and we 
need to present the observed velocities in this manner. 



2.2.2 Long-slit spectra kinematics 

Part of our kinem atical data are obtain ed by means of long- 
slit spectroscopy IjCoccato et al.l 120091 ) from spectra taken 
along the major and the minor axes of the galaxy. The 
width of the slit was 1". Thus we have along each axis pro- 
files of the rotation, of the velocity dispersion, as well as 
two Gauss-Hermite moments (/13 and /14, see appendix |A")) . 
When constructing the models we will use the rotational ve- 
locity profiles along the major axis and of the velocity dis- 
persion along both principal axes (while the Gauss-Hermite 
moments will be used for analysis of the so constructed mod- 
els, see Sect. [4]). We use these data almost "as is" without 
any parametrization. We only bin the data suitably. 

Let us describe how we prepare t he profile of th e mean 
velocity along the major axis. From ICoccato et al.l (j2009i ) 
we get the table containing values of the mean velocity for 
different points along the major axis. So we have a set of 
pairs Xi, Vi, where Xi is the coordinate of the data point 
on the major axis and Vi is the observed line-of-sight mean 
velocity at this point. At first we "reduce" these data to 
the first quadrant (see previous Sect.). For the mean veloc- 
ity this implies multiplying Xi and Vi by -1 for each data 
point with Xi < 0. We need to convert these data into a set 
of two dimensional areas with known mean velocity. These 



2 These are only the major and minor axes of the projected im- 
age, and the real principle axes of the galaxy can, of course, be 
different. 



long-slit data were obtained from rather narrow zones along 
the major axis. We assign these data to a wider zone. We 
assign the data along the major axis to the area defined by 
x < 0.5" || if < 10°, where ip is the angle between the cur- 
rent radius vector and the a;-axis measured on the XZ plane, 
and 1 1 is a logical OR. We divide this area along the :r-axis 
into the pieces so that each long-slit data point corresponds 
to a single piece. We do this as follows. We sort the data 
points by Xi. For each data point we define two values as 
U = (xi~\ +Xi)/2 and ut = (xi + x i+ i)/2. We set h = and 
u n = 2x„ — u n -i, where n is the number of data points. For 
each data point we assign a two-dimensional area defined 
as x > h&i&ix < Uik.k{x < 0.5"\\ip < 10°), where && 
is a logical AND. We then join some of these areas in the 
following way. We create an iV-body system with a surface 
distribution of particles obtained from the surface photom- 
etry (see Sect. l2.1[l . The number of particles is N — 500 000. 
For each area we calculate the number of particles which are 
situated in this area. If the number of particles in some area 
is less than 1000 then this area is joined with a neighbouring 
area. The value of the velocity in each composite area is cal- 
culated as the mean value of the velocities of the constituent 
areas weighed by the number of particles in each area. This 
way of binning the observational data is rather unusual but 
it is convenient for us, because we make sure that in each 
area there is a sufficient number of particles. 

The other profiles are prepared in the same way. The 
data along the minor axis are assigned to the area defined 
as z <0.5"||<£ >80°. 

In our algorithm we don't use any information on the 
errors of the data. But we will use binned long-slit data for 
figures, so, in order to plot error bars, we need to calculate 
the errors of the binned data. We calculate the error of a 
binned datum as e = — i^ X] e *i where n is the number of 
original data points from which binned datum was calcu- 
lated, and e; are the errors of these data points. 



2.2.3 PNe kinematics 

We use the PNe data in the outer part of the galaxy, which 
long-slit data can not reach. More precisely, we use all PNe 
whose distance from the galaxy center R xz — \/x 2 + z 2 is 
larger than 124" (9.5 kpc). 

We "reduce" the PNe data to the first quadrant (see 
Sect. 12.21) . We divide the first quadrant into five zones 
with an equal opening angle ip and use the p angle of the 
individual PNe to place them into the appropriate zone. 
The parameters of these five PNe groups are shown in ta- 
ble [l] For each group we define an area ((p' mln ' < ip < 
<^ (max) ) && {Rif n) < R*z < i?i? ax) ). These areas are shown 
schematically in Fig. [T] and we will refer to these areas as 
Al, A2, A3, A4, A5. 

For each group (area) we calculate the mean velocity 
and the velocity dispersion with the corresponding errors. 
We note that here we take the standard deviation as the 
dispersion. For the calculatio n of the standard deviation we 
use an unbiased estimator (IKennev and Keeping Il95ll . p. 
171), although for such relatively large samples (n > 10, see 
table [1} , the use of an unbiased estimator is not essential. 

We want to draw attention to two features of the PNe 
velocity distribution. The first feature is that the velocity 
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Table 1. Parameters of the five planetary nebulae groups. The first column gives the name of the corresponding area (see Fig.[TJ and the 
second and third ones give its lower and the upper azimuthal boundaries (</j( mm ) and (^( max ), respectively). The remaining columns give 
parameters relevant to the group of PNe in the corresponding area. Here n is the number of PNe, Ri^' n and i?£ are the minimal 
and the maximal distance to the galaxy center, V is the mean velocity, AV is the error of the mean velocity, a is the velocity dispersion 
(standard deviation) and Act is error of the velocity dispersion. 
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3 METHOD 



3.1 General outline 
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Figure 1. Five PNe areas. All planetary nebulae have been 
moved into the first quadrant. 



dispersion at intermediate angles is noticeably higher than 
the velocity dispersion along the major or the minor axes. 
The velocity dispersion in area A3 is twice as high as in ar- 
eas Al or A5. This feature can be clearly seen on the right 
panel of figure 7 in N09, where, except for the inner areas 
the iso-dispersion contours are elongated in a direction in- 
termediate between the two principal axes. Since this figure 
was constructed under the assumption of triaxial symmetry, 
this peculiarity remains true independent of the axisymmet- 
ric assumption. The second feature is that the value of the 
mean velocity in A5 is negative and that it is low in the inter- 
mediate areas A2, A3, and A4. This feature can be explained 
by a twisting of the rotation axis in the outer part of the 
galaxy (see left panel on figure 7 in N09 and corresponding 
discussion), but can not be reproduced by the models with 
an axisymmetric velocity distribution. We, therefore, can 
not take into account in the iterative process the mean ve- 
locity from all regions. Since in axisymmetric systems there 
is no rotation along the minor axis and since the value from 
the Al area agrees very well with the long-slit data (figure 
3.d), we will choose to use the mean velocity from area Al. 



We want to construct an equilibrium A r -body model of an 
elliptical galaxy from its observational data. As described 
in the previous section, for the galaxy under consideration 
we have surface photometry, a distance estimate and various 
line-of-sight kinematics. Assuming that the stellar mass-to- 
light ratio (M/L) is constant, we can obtain from the surface 
photometry and the distance the surface mass distribution 
to within an unknown multiplicative constant M/L. In the 
case of an A^-body model, this implies that we have the pro- 
jected surface distribution of particles, but the mass of the 
individual particles is not known (in our models all particles 
have the same mass). We note, however, that the M/L can 
not be arbitrary, because it is related to the line-of-sight ve- 
locity dispersion in the central part of the galaxy which we 
know from observations (see Sect. 13.4) . We also assume that 
the galaxy is axisymmetric. Observations do not give us the 
inclination of the rotation axis, so this is a free parameter. 

As result, our task is to construct an equilibrium N- 
body model with the given projected surface distribution 
of particles and the given line-of-sight kinematics. The total 
mass of the model is unknown (M/L is unknown), but should 
be found somehow. And the model should be axisymmetric 
with a given axis of rotation. The last condition is of course 
optional, but simplifies the modeling. 

In RAS09 we presented an iterative method for con- 
structing equilibrium A^-body models with given properties. 
The idea of the iterative method is very simple. It relies on 
constrained, or guided, evolution. We simply evolve the sys- 
tem while constraining the desired system properties (see 
RAS09). The same conceptually method, with only rela- 
tively minor modifications can be applied to our present 
task. 

The scheme of the modified iterative method which we 
use in this work is outlined schematically in Fig. [2] We will 
first briefly describe the whole method and then each part 
of the algorithm in detail. We start by building an initial 
A^-body model which has a rigid halo and the desired pro- 
jected surface distribution of the particles. The distribution 
of particles along the line-of-sight can be arbitrary. The ve- 
locities of the particles and the total mass of the system 
(M/L) can also be arbitrary. This model will be the starting 
point for the iterative procedure. At the start of each iter- 
ation we calculate the evolution of the system over a short 
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Figure 2. The scheme of the iterative method which we use for 
constructing the TV-body model of an elliptical galaxy from ob- 
served data. 



detail in RAS09 (section 2.4). We note that this algorithm 
has a free parameter n n b — the number of neighbours. 

It is very easy to modify this algorithm to fix only the 
projected surface distribution of particles. We need to con- 
struct a "new" iV-body system with a given surface dis- 
tribution of particles. In the new model only the x and z 
coordinates of the particles are defined. The y-coordinates, 
the velocities and the mass of the particles should be carried 
over from the old model. The total mass of the new system 
is set equal to the total mass of the old system. We note 
again that in our models all particles have the same mass. 
In the algorithm described in RAS09 we "transfer" from the 
old to the new model only the velocities. In the present case 
we need to "transfer" also the y-coordinates of the particles. 
When we search for the nearest particle we need to do so in 
the two-dimensional space X — Z. The y-coordinate should 
not be taken into account because we copy it from the old 
model particle to the new model particle together with the 
velocities. 



time. Then we put the model through a specific procedure 
which adjusts the total mass of the system (see Sect. 13. 4[) 
every second iteration. We then, if necessary, impose the 
condition of axisymmetry about the given axis. We do this 
in the usual way, i.e. by randomizing the particle azimuthal 
angles. Next, we fix the surface distribution of the particles 
(see Sect. 13.2]) . Finally, we fix the line-of-sight kinematics to 
agree with the observations (see Sect. [33)) . We repeat this 
iteration cycle a number of times until the velocity distri- 
bution and the total mass of the system stop changing. The 
model at this stage is already in equilibrium, or very close to 
it. To be sure of this, we calculate the evolution of the sys- 
tem over a long time scale (3 Gyr in this article) after which 
we obtain the final model in stable equilibrium. Let us stress 
that we consider this long time scale evolution only as a part 
of our algorithm for constructing equilibrium models. 

In the axisymmetric case we assume that the galaxy 
rotates about the axis which lies in the Z — Y plane. This 
rotation axis is defined by its angle a with the Z axis. If 
a = then the rotation axis coincides with the Z axis, so 
the galaxy is edge-on. 

Having described the general outline of our method, we 
will now describe each individual step. 

3.2 Fixing the projected surface distribution of 
particles 

Here we will describe how we fix the projected surface distri- 
bution of particles. We do this using a method very similar 
to that described in RAS09 for fixing the particle mass dis- 
tribution. The idea of that algorithm is as follows: We start 
with a iV-body system, which is the result of a short evo- 
lution (Fig. and to which we will refer to as the "old" 
system. We need to fix the mass distribution in this model 
according to given constraints. We create a "new" TV-body 
system with the desired mass distribution, and we "trans- 
fer" the velocity distribution from the "old" to the "new" 
model. The basic idea of our velocity transfer algorithm is 
very simple. We assign to the new-model particles the veloc- 
ities of those particles from the old model that are nearest 
to the ones in the new model. This algorithm is described in 



3.3 Fixing the line-of-sight kinematics 

Let us first describe an algorithm for fixing the line-of-sight 
mean velocity for the case where we do not assume any sym- 
metry in the system. Our task is as follows. We have an 
TV-body system and some two dimensional area on the sky 
plane where we need to fix the line-of-sight mean velocity 
to the observed value. In our models the sky plane coincides 
with the XZ plane, so the line-of-sight velocity is the ve- 
locity along the Y-axis. The given area is a two dimensional 
area in the XZ plane, so when we search the particles which 
belong to the given area we do it regardless of y-coordinates 
of the particles. We denote by v y the desired value of the 
line-of-sight mean velocity in the given area and by v' y the 
mean value of the y velocities of all particles in the area. We 
need to change somewhat the particle velocities so that v' y 
becomes equal to v y . This is achieved by setting the new y 
velocity component of particle i in the given area to 

Vyi = V U + (Vv ~ K) ' 



where v' yi is the current value of the i- 



(2) 

th particle and v y i is 
the corrected i-th particle y velocity. 

The algorithm for fixing the line-of-sight velocity disper- 
sion is very similar. Let us denote by a y the desired value of 
the line-of-sight velocity dispersion in the area under consid- 
eration. In the given area we calculate the current value of 
the line-of-sight velocity dispersion a' y and the current value 
of the line-of-sight mean velocity v' y . The new y velocity 
component of particle i in the given area is set as 



+ v v 



(3) 



Let us now consider a galaxy with the following symme- 
tries (as described in Sect. 12. 2) . The galaxy is axisymmetric, 
with a symmetry axis which lies in the YZ plane, and also 
has a reflection symmetry with a plane of symmetry perpen- 
dicular to the axis of symmetry and containing the center of 
coordinates (0, 0, 0). In this case, if we know the line-of-sight 
velocity distribution in the single quadrant then we know it 
for the whole system (see Sect. 12.21 for details). So, if we as- 
sume such a symmetry then all observed kinematical data 
can be "reduced" to the first quadrant (x > 0, z > 0). 
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In the case of such a symmetry the task of fixing the 
line-of-sight kinematics is as follows. We have some two di- 
mensional area in the first quadrant of the sky plane for 
which we have the line-of-sight mean velocity and/or the 
line-of-sight velocity dispersion. We need to fix these kine- 
matic parameters in the given area taking into account the 
discussed symmetries. We invert the sign of v y for each par- 
ticle with x < and then we flip all particles to the first 
quadrant, i.e. we set x = \x\ and z — \z\ for each particles. 
Now we can apply the algorithms for fixing the line-of-sight 
kinematic parameters which we described in the first part of 
this section. We then flip all particles back to their original 
positions and we invert the sign of v y for each particle with 
x < 0. 



3.4 Adjusting the total stellar mass 

Since we have no a priori knowledge of the total stellar mass 
in the galaxy, we have to construct our models assuming 
that the stellar M/L is unknown. In our case this implies 
that the total mass of the system is unknown. This of course 
raises questions. Does an equilibrium model with proper- 
ties in agreement with observations exist for a unique value 
of the mass? Or for a range of values? For example, for 
a given equilibrium, spherical, isotropic model with known 
projected surface distribution of particles and line-of-sight 
velocity dispersion in the center of the model, the total mass 
of the model would be uniquely determined. But in our case 
it is not so obvious. Moreover, as we will show, models with 
different inclinations of the rotation axis have slightly differ- 
ent total masses (see Table [2 for example models ALO and 
AL45). So if we don't fix the inclination of the rotation axis 
then the total mass may well not be unique. 

Nevertheless we need to find some value of the total 
mass for which we can construct the equilibrium model. The 
straightforward way to solve this problem is to construct a 
series of models with different total masses and choose the 
one which, according to some definition, is closest to equi- 
librium, or to find ranges of values for which the resulting 
equilibrium model is in agreement with the observational 
constraints. This, however, would be excessively time con- 
suming, due to the large number of models that need to be 
constructed. 

We will, therefore, use a different algorithm, in which we 
adjust the total mass during the iterative process (Fig. [2]). 
As we described earlier, in each iterative step we let the 
system evolve on a short time scale (Fig. [2]). We calculate in 
the beginning and in the end of this evolution the velocity 
dispersion along the line-of-sight in some given part of the 
system, which we denote by ci and 02, respectively In all 
experiments described here this given part was a sphere with 
radius equal to 10 kpc. We note that the value of o\ is partly 
defined by the given line-of-sight velocity dispersion which 
we fix on the previous stage of the iteration (see Fig.[2|. On 
the other hand the value of <72 is influenced by the total mass 
of the model. So, if we choose a value of total mass which is 
not appropriate and we don't change it during the iterative 
procedure, then ci and 02 will be different after any number 
of iterations. We want to construct the equilibrium model so 
that the values of g\ and 02 are close, so we adjust the total 
mass in the system by multiplying all velocities by a factor of 
<Ti/(T2 and the masses of all the particles by a factor oia\ja\. 



By trial and error we found that the iterations converge 
faster if we adjust the total mass not at every iteration but 
only every second iteration. 

For a model without dark halo this means that, if the 
system at the end of the evolution was in equilibrium, then 
the rescaled system will also be in equilibrium, but will have 
the line-of-sight velocity dispersion in the selected part of the 
model equal to o\. For models with halo this explanation 
is of course not valid, but this is not a serious drawback 
since there will be further iterations to bring the system to 
equilibrium. We thus used this algorithm for constructing all 
our models, including models with dark halo, and it worked 
well in all cases. I.e. all these models constructed by means 
of the iterative method with this mass-adjusting algorithm 
were in equilibrium. But we note that it is possible that, 
for models with a very massive halo dominating the central 
part of the galaxy, this algorithm may not work and the 
iterations would not converge. In such cases one would have 
to resort to the straightforward but very time consuming 
algorithm described above. 



4 MODELS 

4.1 Description of the models 

We constructed models with three types of halo. The first 
type is models wit hout halo. The second and third type are 
rigid NFW halos (|Navarro et al.l Il99rj . I1997T ) with density 
profile 



p(r) = 



Ps 



(r/r s )(l + r/r s ) 2 



(4) 



where p s and r s are the characteristic density and scale ra- 
dius of the halo. The second type was found in N09 to be 
the best-fitting NFW halo for this galaxy (see Sect. 4.2.5 
in N09). This is a relatively light halo with parameters 
p a — 0.0019 Mqpc~ 2 and r s — 32 kpc, i.e. a concentration 
parameter c vlI ~ 8 and a virial mass M v i T « 10 12 Mq (see 
N09) . The third type is a relatively massive halo with param- 
eters p a = 0.00522 M e pc~ 2 and r a = 26.5 kpc. N09 found 
that a model with such massive halo is a relatively poor fit 
of the observational data (see their Sect. 4.2.5). Since we use 
a totally different approach, we include this model to test 
whether this conclusion is method-dependent, or not. This 
halo has c vir = 12.3 and M vir = 2 x 1O 12 M . 

Another parameter of our models is the angle a defin- 
ing the inclination of the rotation axis (see Sect. [3}. We 
construct models for two values of the angle a, namely 0° 
and 45° . 

As we discussed in Sect. l2.2T3l the azimuthal variation of 
the PNe velocity dispersion presents an interesting feature 
in the outer parts. Namely, the velocity dispersions in areas 
close to the principle axes (areas Al and A5) are consider- 
able smaller than in intermediate areas (areas A2, A3 and 
A5). It is not clear whether such a feature can be reproduced 
by an equilibrium axisymmetric system. So we construct two 
sets of models. In the first ones, which we denote as "A" , we 
don't try to model this feature. Thus, we ask the iteration 
method to fit the following quantities: 

• The projected surface distribution of the particles (see 
Sect. Ell 
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• The radial profile of the mean velocity along the major 
axis obtained from long-slit spectroscopy 

• The radial profiles of the velocity dispersion along the 
major and minor axes obtained from long-slit spectroscopy 
(see Sect. l2^2jl 

• The mean velocity in area Al obtained from PNe kine- 
matics. Note that we don't fix the mean velocity in other 
areas (see Sect. 12. 23]! 

• The velocity dispersion in the two PNe areas Al and 
A2, which are close to the major and minor axes, respec- 
tively. 

For the second set of models, which we denote by "B", 
we take into account all available information. We, there- 
fore, try to fit the velocity dispersion in all five PNe areas 
separately. This introduces a considerable extra difficulty be- 
cause of the peculiar PNe velocity distribution. As already 
mentioned, such a fit may of course not be possible with 
models such as ours, i.e. with models that are both axisym- 
metric and in equilibrium, ft is, nevertheless, important to 
try for at least two reasons. First it is useful in all cases to try 
and fit all available data, since only that can tell us how far 
off the attempted fit is from reality. Second, since we allow 
the model total freedom regarding the velocity anisotropy, 
a fit may be found, in which case an interesting effect of 
anisotropy would be revealed. 

We will, therefore, discuss in total 12 models, which we 
will denote as follows. The first symbol in the name of the 
model denotes the set (A or B) and the second denotes the 
halo. Here "N" is for models without halo, "L" is for models 
with a relatively light NFW halo (p s = 0.0019 M Q pc' 2 , 
r s = 32 kpc) and "H" is for models with a relatively massive 
NFW halo {p s = 0.00522 M Q pc~ 2 , r s = 26.5 kpc). The 
last number in the model name denotes the angle a of the 
rotation axis inclination. 

Let us now discuss in some detail how we construct 
these models. There is considerable freedom in choosing the 
initial model, from which the iterative search will start. Our 
basic initial model is a model with the given surface distri- 
bution of particles (generated using the rejection method), 
and with zero velocities. The y coordinates of the particles 
were chosen as random numbers from the interval [—1,1] and 
the initial mass was chosen according to M/Lv — 1. These 
choices are of course totally ad hoc, but this does not matter 
since other values, although starting the iterative procedure 
(see Fig. [5| from a different initial model, lead to essen- 
tially the same final model. In practice, there is only one 
significant property of this initial model, namely whether it 
is rotating (or not). Indeed, if we start the iterative search 
from a rotating, rather than a non-rotating initial model, 
we will end up with a model that is different, albeit not in 
all properties. A rotating initial model can be constructed 
as follows. We first take the basic non-rotating initial model 
and put it in the iterative procedure (see Fig. [2]| . After a rel- 
atively small number of iterations, the model becomes close 
to steady state. We then choose the axis of rotation and 
set, for all particles, the azimuthal velocity with respect to 
the chosen axis equal to the circular velocity. We have thus 
obtained a rotating initial model. 

Models with a — 45 constructed from rotating and from 
non-rotating initial models are practically identical. Non- 
inclined models (a = 0), however, are slightly different. In 



Table 2. Stellar M/L and relative halo mass for the constructed 
models. The first column gives the name of the model, the sec- 
ond and third columns the stellar mass-to-light ratios in V and 
B bands, M/Ly and M/Lb, respectively. The fourth and fifth 
columns give the ratio of the halo mass, M^ix), to the stellar 
mass, M,(i), both calculated within a sphere of radius x = R e 
and x = 5iJ e , respectively. 



Model 


M/L v 


M/L B 


M h (R e ) 
M,(R e ) 


M h (5R e ) 
M„(5R e ) 


AN0 


4.16 


4.82 








ALO 


3.55 


4.11 


0.13 


0.96 


AH0 


2.76 


3.22 


0.36 


2.52 


AN45 


4.23 


4.90 








AL45 


3.73 


4.32 


0.12 


0.90 


AH45 


3.15 


3.65 


0.30 


2.18 


BN0 


4.21 


4.89 








BL0 


3.55 


4.12 


0.13 


0.96 


BH0 


2.89 


3.35 


0.35 


2.41 


BN45 


4.27 


4.95 








BL45 


3.81 


4.41 


0.11 


0.88 


BH45 


3.19 


3.69 


0.30 


2.16 



particular, there is a difference in rotation in areas far from 
major axis. We note that we fix the mean line-of-sight ve- 
locity only in areas close to major axis (long-slit data along 
major axis and mean velocity in Al area). It is not surpris- 
ing that models constructed from initially rotating models 
rotate slightly faster in areas far from major axis. But this 
difference is not sufficiently significant to warrant further 
discussion. 

Here we will discuss models constructed from the ro- 
tating initial model. Each model was constructed as follows. 
At first we construct a basic non-rotating initial model with 
N = 300 000 particles. We put this model into the iterative 
procedure and make 100 iterations with relatively low pre- 
cision. This is possible, because each iteration is very short 
and errors do not acc umulate (RAS09). W e use the fast N- 
body code gyrfalcON (|Dehnenll2000al . 120021 ') with an integra- 
tion step and a softening length equal to dt — 1/2 12 Gyr and 
e = 0.05 kpc, respectively. The tolerance parameter for gyr- 
falcON was set to 6t — 0.9 and the duration of each iteration 
to ti — 0.05 Gyr. Using this constructed model, we create 
the rotating initial model with N = 500 000 particles. Again 
we put this model into the iterative procedure and make 
500 iterations. The integration step and the softening length 
were taken dt — 1/2 14 Gyr and e = 0.02 kpc, respectively. 
The duration of each iteration and the tolerance parameter 
for gyrfalcON were chosen as in the previous stage. The fi- 
nal stage of our procedure is a free evolution over a time 
scale of 3 Gyr (see Sect[2]), with parameters dt = 1/2 Gyr, 
e = 0.02 kpc and 9± = 0.6. We would also like to mention 
that before fixing the parameters to these values, we made a 
number of tests, such as increasing the number of particles 
up to tenfold, and did not find any significant improvements 
in the fits. 



3 We use a system of units where the unit of length is u; = 1 kpc, 
the unit of velocity is u v = 1 km/sec, the unit of mass is u m = 
10 10 Mq and consequently the unit of time is ut ~ 0.98 Gyr. For 
simplicity, all time values in this papers are presented with the 
assumption that ut = 1 Gyr. 
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Figure 4. Radial profiles of velocity anisotropy for models AN0, 
ALO and AHO. The velocity anisotropy parameter is calculated 
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as p = 1 §■ , where a r and cr^, are the velocity dispersion in the 

radial and the <p direction, respectively, in a spherical coordinate 

system. 



The stellar mass-to-light ratio and relative mass of the 
dark halo inside one and five effective radii, respectively, are 
shown in Table[2] According to N09 the value of the effective 
radius R e — 48". 2 ~ 3.69 kpc. We use photometry in the V- 
band, so for us it is more straightforward to calculate the 
mass-to- light fraction in V-band M/Lv- To compare our 
results with N09 we also calculate mass-to-light fraction in 
B-band M/ Lb- We assume that (B — V) is 0.65 for the Sun 
and 0.81 for NGC 4494 (see N09). We note that, if we change 
the adopted distance, the mass-to-light will be changed (see 
Sect. [2]). For all models, there is considerably less halo mass 
than stellar mass within lR e . For models with a "light" halo, 
the mass of the dark matter inside 5R e is approximately 
equal to the mass of the stellar component, while for model 
with a "heavy" halo the mass of dark matter inside 5R e is 
more than twice that of the stellar component. 



4.2 Discussion of the first set of models 

4-2.1 Basic comparison 

When constructing our first set of models, models A, we use 
only velocity information from the vicinity of the principle 
axes. In particular, we take into account the velocity dis- 
persion in PNe areas Al and A5 and ignore observations 
in areas A2, A3, A4. In so doing we want to demonstrate 
that our method can be used to construct galactic models 
following given observational constraints, both photometric 
and kinematic. 

Let us first discuss models AN0, ALO and AHO, which 
have an inclination angle a — 0, i.e. their rotation axis is 
parallel to the sky plane. Our results for these models are 
summarized in Table 2 and Figs. [3] and [4] In order to increase 
the resolution and reduce the noise in our figures we use 



a trick described by lAthanassoulal (|2005bj). We stack ten 
snapshots closely spaced in time, with At = lOMyr, and we 
calculate all values for this combined snapshot. This allows 
us to reduce the noise in our plots quite significantly. 

Since we include rotation, we compare models to obser- 
vations separately for mean velocities and for velocity dis- 
persions, and do not fold these two quantities into a single 
root-mean-square velocity v rms = y/v 2 + o~ 2 (N09). 

In general, there is a good agreement between the mod- 
els and the observational data which were used by the it- 
erative method for constructing them (Fig. [3}. We get an 
excellent fit of the density profile (panels (a) and (b)). We 
also get good agreement with the mean velocities on the ma- 
jor axis except for radii within the centermost region (panels 
(d) and (e)) and with the velocity dispersion on major and 
minor axes (except for the bump at radii between 60" and 
90" - see panels (f) and (g)). All these agreements hold for 
all three models, i.e. models AN0, ALO and AHO. Also all our 
model fit very well mean velocity of PNe in area Al (panel 
(h)). Velocity dispersion in PNe areas will be discussed later 
in this section. Let us now discuss in more detail the parts 
where the models fail to reproduce the observations. 

• As many other ellipticals, NGC 4494 has a kinemati- 
cally decoupled core. This is clear not only from the mean 
velocity, but also from the velocity dispersion profiles along 
the major and minor axes. Such features are believed to be 
out of equilibrium, e.g. due to a merger with a a small com- 
panion that reached the galactic centre by dynamical fric- 
tion. Furthermore, the material in this region may not be ax- 
isymmetrically distributed and/or have a different orienta- 
tion from the remaining galaxy. Arguments in favour of these 
possibilities are that the density distribution near the cen- 
ter does not follow the Sersic law and that the profile of the 
mean velocity does not have the shape expected for axisym- 
metric mass distributions in equilibrium. Departures from 
axisymmetry and from equilibrium are beyond the scope 
of this paper. It is thus normal that neither N09 - with a 
spherically symmetric non-rotating equilibrium model - nor 
we - with an axisymmetric, rotating or non-rotating equilib- 
rium model - reproduce the structure in the innermost part. 
Models including kinematically detached cores are beyond 
the scope of both studies. 

• The velocity dispersion profiles have a bump along both 
the major and minor axes in the area between 60" and 90" 
(panels (f) and (g) of Fig. [3|. N09 considered the long slit 
data as kinematical constraints only up to 60", where they 
considered them as more accurate. But their plots show 
clearly that their models have no such bump in the 60" 
to 90" region. We kept all the long-slit data as constraints, 
but none of our models reproduced this bump. It is possible 
that in NGC 4494 this bump is transient, presumably due 
to some collision event and thus can not be reproduced by 
the Jeans method, or by our iterative method, which can 
only find equilibrium solutions. This, together with the pre- 
vious discussion on the centre-most area, argues that some 
regions of NGC 4494 are not in exact steady state, as often 
observed and as could be expected in the framework of the 
ACDM paradigm. 

• The ellipticity of model AN0 agrees on average with the 
mean observed value (panel (c)). However models ALO and 
AHO have bigger ellipticity values than observed in their 
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Figure 3. Comparison of models ANO, ALO and AHO with the observational data, (a) and (b) show the dependence of n xz on R xz , 
where n xz is the number of particles in concentric cylindrical shells and R xz = yjx 1 + z' 1 . The thick solid line shows the profile calculated 
for a model with the input data prepared as described in Sect. 12.11 Panel (c) shows the ellipticity profiles of the models projected on 
the sky plane (XZ plane), calculated with the IRAF task ELLIPSE. Panels (d) and (e) show the profile of the mean velocity along the 
major axis and panels (f) and (g) the profiles of the velocity dispersion along the major and minor axes, respectively. In (d) and (f), 
together with the long-slit data, we show PNe data in area Al which is close to the major axis (see Fig. [TJ. The value calculated for this 



area is assigned to radius (R 



(mm ) 



R 



(max) 



)/2 (Table [!}• In panel (g), together with the long-slit data, we show the PNe data in area 



A5 which is close to the minor axis. In panels (d), (f) and (g) the observed data are binned as described in Sect. l2~2.2l Panel (h) shows 
the mean line-of-sight velocity calculated for the five PNe areas (see Sect. l2~2.3H . Note that only the mean velocity in area Al is used for 
constructing the models we use. Panel (i) shows the line-of-sight velocity dispersion calculated for the five PNe areas. Only dispersions 
in areas Al and A5 are used for construction of A models. The radial profiles for models ANO, ALO and AHO are given by solid, dashed 
and dotted lines, respectively. The open squares show the long-slit data, not binned for panel (e), and binned as discussed in Sect. ?? 
for the remaining panels. The filled circles show the PNe data, which, following N09, were used as observational constraints for these 
models, and the X signs the remaining PNe data. 
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Figure 5. As in Fig. [3] but for models AN45, AL45 and AH45. 



outermost parts. This is partly due to the fact that these 
two models have slightly boxy isophotes at their periphery, 
which we believe to be connected to the very high veloc- 
ity anisotropy in their periphery (Fig. [4J, i-e. to be linked 
to the radial orbit instability. As already mentioned, our 
algorithm for model construction includes a self-consistent 
A^-body evolution on a relatively long time scale (Fig. [2|, so 
that our models can "feel" such instabilities. Let us, how- 
ever, underline that this effect is relatively small and the 
isophote boxyness is rather subtle. As a result the difference 
with the average ellipticity is of the order of only 0.05, i.e. 



much smaller than the corresponding difference for spherical 
models, such as those of N09. 

We can conclude that our models reproduce well the 
projected surface density, the mean apparent ellipticity, the 
mean velocities and their dispersions (the latter though not 
over the full radial and angular extent of the galaxy). The 
corresponding fits are no worse than the N09 models. The 
discrepancies between our models and observations occur 
mainly in regions whose data were not considered by N09. 

Fig. 0] shows the radial profile of the velocity anisotropy 
for our a — models and clearly illustrates the degeneracy 
between mass and anisotropy. These three models have dif- 
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ferent halo masses, but the same stellar density profile and 
observed kinematics (Fig. [3} and this is possible because 
they have different velocity anisotropics. As was expected, 
radial anisotropy is increasing with halo mass (see Fig. [4} . 

The inclined models AN45, AL45 and AH45 are con- 
siderably different from the previously discussed ANO, ALO 
and AHO models, as they have an intrinsic ellipticity of 0.35. 
The quality of the fits, however, is very similar (Fig. [5}. 

Seen the size of the error bars and the small effect of 
the halo mass on the density and velocity radial profiles, 
it is not possible from the above comparisons alone to rule 
out any of the models, or even to set a strong preference 
to one rather than another. The only exception is the el- 
lipticity profiles, which are considerably better fitted by the 
halo-less models ANO and AN45 (Figs. |3]and[5]). On the 
contrary, all models with dark halo have in their outer parts 
slightly boxy isophotes (presumably due to their higher ve- 
locity anisotropy) and a higher than mean observed elliptic- 
ity. In our modeling, however, we use only the mean value 
of the ellipticity and not the full profile. Also we don't use 
high order isophote shape parameters, so we enforce our 
models to have precise elliptical isophotes. NGC4494 has a 
negative 0,4 isophote shape parameter on periphery, i.e. boxy 
isophotes (see appendix A in N09). So we cannot rule out 
models with dark halos relying only on the fact that they 
have slightly boxy isophotes in their outer parts. 

Models without dark halo fit perfectly the velocity dis- 
persion of the PNe that have been used as observational 
constraints, i.e. in regions close to the major and minor axes 
(see panels (i) of Figs. [3] and |5J. Model AL45, however, with 
a relatively light dark halo, also fits the observational data 
along the major and minor axes rather well, while ALO fits 
only slightly worst (see panels (i) of Fig. [3] and [SJ). And, seen 
the error bars, even model AH45 with a relatively massive 
dark halo can not be ruled out. 



4-2.2 Higher order moments of the velocity distribution 

When we constructed the models, we did not use the higher- 
order moments of the velocity distribution as observational 
constraints. It is thus interesting to compare the third and 
fourth order moments of the velocity distribution of our 
models to those of the observational data. Moreover, such 
moments may help break the mass-anisotropy degeneracy 
and were used by N09 to argue for the need of a low mass 
halo. At the periphery of the galaxy we have kinematical 
information only from the PNe. We thus calculate the skew- 
ness and the kurtosis in the five PNe areas (Fig. [TJ and com- 
pare them to the corresponding model values, calculated for 
our models as described in appendix [B] It should, however, 
be noted that the number of planetary nebulae is very small 
and that this affects more the higher order moments. Thus, 
even assuming a Gaussian velocity distribution, the uncer- 
tainties for skewness and kurtosis are very big, while without 
this assumption the uncertainties are formally infinite (ap- 
pendix [B}. 

One of the advantages of our models is that we can 
easily calculate any parameters for them. This will allow us 
now to calculate for our models exactly the same high-order 
moments of the velocity distribution as for observations. We 
can thus calculate two Gauss-Hermite moments (/13 and hi) 
along the major and minor axes of the galaxy, as described 



in appendix and compare them with the corresponding 
profiles from long-slit spectroscopy. The results are given in 
Figs. IS] and [7J for a — and a — 45, respectively. These fig- 
ures show clearly which of the higher order moments depend 
on the halo mass and in which way. 

The most interesting of the Gauss-Hermit moment pro- 
files is probably the /13 profile along the major axis. This 
is visibly different for models with different halo mass (see 
panels (a) of Figs. [5] and [7]), the most massive halo one hav- 
ing the highest /13 values and ANO the lowest. The best fit 
seems to be for the intermediate halo mass. To establish this 
we calculated the x 2 > excluding the region of the kinemat- 
ically decoupled core, and normalised it by the number of 
data points. We find that for a = 0, the values are 3.9, 0.4 
and 5.7 for models with no halo, light halo and heavy halo, 
respectively. The corresponding numbers for a = 45 are 2.4, 
0.8 and 1.8, respectively. These numbers show a preference 
for the models with light halo, and argue, albeit weakly, for 
a preference for a — 0. The profile of /14 along the major axis 
also shows some dependence with halo mass, but much less 
so than the corresponding /13 profile, (panels (d) of Figs. [5] 
and[7J). Profiles along the minor axis for both Gauss-Hermite 
moments show no clear dependence on halo mass (panels (b) 
and (e) of Figs. |S]and[7|, except for the hi profiles in the 
outermost parts, where, however, there are no long-slit data. 
In the central part of the model, i.e. inside ~ 30", all our 
models have /14 values which are clearly less than those of 
the observations for both axes (Panels (a) of Figs. [6]and[7]l 
and this is true also for the /13 major axis profiles. This could 
again be linked to the kinematically decoupled core. 

The kurtosis values calculated in each of the PNe areas 
for all "A" models are, for practical purposes, almost the 
same (panels (f) in Figs. [S] and [7J, the differences being 
much smaller than the observational errors. The situation 
with the skewness is the same (panels (c) in Figs. [6] and 
[JJ. We, therefore, do not believe that the comparisons in 
panels (c) and (f) of Figs.|6]and[7]can be used to distinguish 
between models of different mass. 

To summarise, the /13 profiles along the major axis pro- 
vide some arguments in favour of models with a light halo 
(particularly ALO and AL45), in good agreement with what 
was found by N09. 

4.3 Discussion of the second set of models 



As discussed in section 12.2.31 the velocity distribution of 
the PNe has an interesting feature, namely that the veloc- 
ity dispersion near the principal axes (regions Al and A5) 
has much smaller values than the velocity dispersions in the 
intermediate areas (regions A2, A3 and A4). When build- 
ing the models presented in the previous section, we did not 
use this as an observational constraint and ignored observa- 
tions in intermediate areas. As results all our "A" models 
have almost the same dispersion in all PNe areas (see pan- 
els (i) in Figs. |3] and [5]l, so they obviously fail to reproduce 
observational data in intermediate areas. Here we consider 
more realistic models ("B" models) for whose construction 
we used all available observational data, to see whether the 
tangential variation of the velocity dispersion can be repro- 
duced by rotating axisymmetric models. 

In general both sets of models are very similar, with 
differences only at the periphery of the galaxy. Furthermore, 
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Figure 6. Comparison of higher moments of the velocity distribution of models AN0, AL0 and AH0 with the observational data. Panels 
(a) and (b) show profiles of the /13 Gauss-Hermite moment along the major and minor axes, respectively (calculated as described in 
appendix [X} . Panel (c) shows an estimate of the skewness Gi calculated for the five PNe areas (calculated as described in appendix fBll , 
(d) and (e) show profiles of the /14 Gauss-Hermite moment along the major and minor axes, respectively, (f ) shows an estimate of the 
kurtosis G2 calculated for five PNe areas. These estimations of the model skewness and the kurtosis are "reduced" to a sample size equal 
to the number of PNe in the corresponding areas (see appendix |Bj . The error bars in panels (c) and (f) were calculated assuming a 
Gaussian distribution function (see appendix |B]| . The symbols and line styles are as in Fig. [3] 



the quality of the fit to the observation data is also similar. 
Thus, to save space, we show only a very restricted set of 
figures. Fig.[8]shows the velocity dispersion calculated in the 
five PNe areas for our B models. 



We find that only non inclined models with a dark halo 
(i.e. models BL0 and BH0) reproduce this feature, and even 
they only partially. Namely, the velocity dispersions near the 
principal axes (areas Al and A5) are too high, while at the 
middle (region A3) it is marginally too low. Thus we were 
able to reproduce qualitatively the form of the azimuthal 
variation, but not quantitatively, since the model amplitude 
is lower than necessary to fit the observations. For these 
models, we can not make specific comparisons with N09, 
because the latter did not try and reproduce this feature. 



We can thus conclude that, although our models re- 
produce most of the observed features well, they fail to fully 
reproduce the velocity distribution of planetary nebulae. We 
will discuss possible explanations in the next section. 



5 CONCLUSIONS 



We used our iterative method 

IjRodionov. Athanassoula fc Sotnikoval I2009T ) to construct 
dynamical models of NGC 4494. One of the advantages of 
our method is that the models are produced in the form 
of TV-body snapshots with equal mass particles, which can 
be directly used in A^-body simulations. This, for example, 
allows us to check directly whether the constructed model 
is in steady state and to calculate any quantities, or 
parameters of these models, which can then be compared 
to observations. As already discussed in Sect. [3] (see also 
Fig- 0) after the main part of the iterative method, we 
let our models evolve freely on a time scale of 3 Gyr, and 
then only do we calculate the quantities to be compared 
with the observational data. Since the aim of this evolution 
is to check whether the models we constructed are indeed 
in steady state, we used the same halo as when building 
the model, i.e. a rigid halo. This evolution showed clearly 
that all our models are in a steady state. For completeness 
sake, however, we also tried simulations with a live halo 
and we find no considerable differences from the rigid halo 
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simulations, arguing that in these cases there is not much 
interaction between the dark and the stellar matter. 

We used the observational data given by N09, namely 
surface photometry, stellar kinematics along the major and 
minor axes as obtained by means of long-slit spectroscopy, 



and velocities and positions of 255 planetary nebulae (see 
N09). We use PNe data only in the outer part of the galaxy 
where long-slit data are absent. The surface photometry 
gives us the surface distribution of particles but not the to- 
tal mass of the system because the mass-to-light ratio is un- 
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known. Our algorithm automatically adjusts the total mass 
of the model (Sect.[3)|. 

We constructed models with three types of halo, all of 
which were already discussed in N09. The first type is models 
without halo. The second type of halo is a relatively light 
NFW halo, which N09 found to be the best-fitting NFW 
for NGC 4494. The third type of halo is a relatively heavy 
NFW halo. We, furthermore, constructed two sets of mod- 
els: models "A" and models "B". These two sets have only 
one difference. When constructing models "A" we used as 
observational constraints only the PNe in areas close to the 
principal axes, neglecting information on the velocity distri- 
bution of the PNe in intermediate areas, on the contrary, in 
order to construct models "B" we used all available infor- 
mation. 

One important goal in making the "A" models was to 
demonstrate the ability of our method to construct equilib- 
rium models with given observation constraints. This was 
fully achieved. In general, there is a good agreement be- 
tween our models and the observational data concerning 
the projected surface density, the mean apparent elliptic- 
ity, the mean velocities and their dispersions (except for 
specific radial ranges, where some non-equilibrium or non- 
axisymmetric structure could be present in NGC 4494). We 
showed that our models reproduce observations data not 
worse than N09 models. But our models have the added ad- 
vantage that they are live, non-spherical, rotating iV-body 
systems. 

A further goal was to see whether it was possible to 
set constraints on the dark halo mass with our models. If 
we don't take into account the high-order moments of the 
velocity distribution, then the best models are the models 
without dark halo. This, however, is only a slight preference 
and models with a light halo and even models with a heavy 
halo can not be rejected. It is thus necessary to try higher 
order moments. We found that the major axis profile of the 
third order moments shows the strongest dependence on the 
halo mass and that the best fit was for the models with the 
light halo. This is in agreement with what N09 found for 
their models, though with a different technique. 

N09 found that a heavy halo gives worst fits to the 
data, without referring to higher order moments, i.e. from 
the root-mean-square velocity profile, which is a combina- 
tion of the mean velocity and the dispersion. This, however, 
is not the case for our models. For example, model AH45, 
which has a heavy halo, fits all low order moments of the 
velocity distribution very well. These two results taken to- 
gether, show that there are no simple spherical models fitting 
the lower order models, but that there can well be more gen- 
eral models that do. Thus, our model which is axisymmet- 
ric, rather than spherical, rotating and inclined with respect 
to the line of sight, has no difficulty with the lower order 
moments. When we place the bar higher, i.e. when we ask 
our models to fit also the third and fourth order moments, 
we find that the models with heavy halos do worse than the 
models with lower mass halos, as we described in Sect. l4.2T2l 
This does not necessarily mean that there are no heavy halo 
models that fit the third and fourth order moments. It can 
simply mean that it is necessary to consider a yet more gen- 
eral model, e.g. a triaxial one, in order to get such a fit. In 
other words, one has to be careful not to extrapolate a given 
result further than the class of models for which it was de- 



rived. And this of course makes it very difficult to set any 
strong constraints to the halo mass in NGC 4494, and to 
ellipticals in general. 

N09 found that some halo is required by comparing the 
kurtosis of their model with that of the PNe velocity distri- 
bution at the outermost parts of NGC 4494. Again, this is 
not the case for our models. It is the hs radial profile along 
the major axis that proved in our case to be more sensitive 
to the halo mass, and thus allowed us to set some preference 
for the light halo model, in agreement with what was found 
by N09 for the fourth order moment. 

When constructing and analysing the "A" models we ig- 
nore the PNe data at angles intermediate between the major 
and the minor axes. Since this is rather ad hoc, we built the 
"B" models, in which all PNe data were used as observa- 
tional constraints. By doing so, we wanted to test whether 
we could make models which reproduced the azimuthal vari- 
ation of the velocity dispersion, i.e. which have relatively low 
values near the principal axes and considerably higher ones 
in between. Although we were able to built models that qual- 
itatively reproduce this feature, we were unable to reproduce 
it quantitatively; i.e our models always display a lower am- 
plitude of this variation than the observations. There are 
several possible explanation to this. 

It should first be noted that the observational data at 
the periphery of NGC 4494 is very sparse. We have only 
around 15 PNe in each of our PNe areas (see Table Q]). So 
even the formal uncertainties of the PNe dispersions in these 
areas are rather big. Moreover, there could still be some 
contamination in the PNe sample (see section 2.2 in N09), 
which would further increase the uncertainties. As a result, 
there is still a considerable probability that the observations 
do not contradict e.g. model BLO (see fig. [3a, error-bars are 
one a ). 

There is of course always the possibility that a model 
which fully reproduces the observations exists, but that 
our method failed to construct it. We made very extensive 
searches and we do not believe that this possibility is likely, 
but, strictly speaking, we can not exclude it. Furthermore, 
the observations point to other, more likely alternatives, 
which we discuss below. 

• In our models the halo is spherical. It is likely that mod- 
els with an axisymmetric or a triaxial halo would reproduce 
observations better. Note, however, that this alternative can- 
not help the models without a dark halo. 

• Our models are axisymmetric, while, as discussed in 
Sect. 12.231 the velocity distribution in the galaxy is not fully 
axisymmetric. It would thus have been better to consider 
more general types of models, for example, triaxial models 
or models with even less symmetry. This, however, would 
introduce further free parameters and is beyond the scope 
of this paper. 

• It is possible that the outer parts of NGC 4494 are not 
in a steady state. As we already discussed, there is evidence 
that NGC 4494 is not exactly in steady state. And if there 
are transient features in the innermost and perhaps in part 
some of intermediate part of the galaxy (the kinematically 
decoupled core and between 60" and 90"), then such fea- 
tures can also be present at its periphery. Moreover, the 
dynamical time scale at the periphery is much longer than 
in the intermediate part, so that any non-equilibrium struc- 
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tures will tend to equilibrium much slower than in regions 
further in. 

• Seen the low number of PNe in the periphery of the 
galaxy, the uncertainies are very large and therefore the con- 
straints in these regions very loose. It might therefore not 
be necessary to dismiss a model if the fits in these regions 
are poor. 

Any of the three first alternatives would necessitate a 
more general model than what we have considered here. 
Triaxial models could, at least in principle, be built, but 
the number of free parameters would increase so that more 
observational constraints would be necessary. However, nei- 
ther our techniques, nor other techniques discussed so far, 
can build non-equilibrium models. This could be done only 
by iV-body simulations, which, however, would have an ex- 
traordinarily unwieldy free parameter space. 

To summarise, we were able to reach our first goal, 
namely we showed that our iterative method (RAS09) can 
be used for building models with given observational con- 
straints. Concerning our second goal, i.e. to set strict con- 
straints on the halo mass of NGC 4494, we can only claim a 
more limited success. Comparing the third order moments 
of our model velocity distribution with the observations, we 
find that it is the light halo model that gives the best fit. 
Even this, however, does not fit all available data, and it 
could be that a more general model (e.g. triaxial) is neces- 
sary. It is, however, not possible to extend the preference of 
light halos to these more general models, without actually 
making them and analysing their properties. It could thus 
well be that it is another type of model (such as with no 
halo, or with a heavy halo) that gives the best fits in such 
a case. Thus the present set of models, although showing 
a preference for a light halo, can not set a very strict con- 
straint to the halo mass of NGC 4494, particularly seen the 
sparseness of the PNe in the outermost parts. 



ACKNOWLEDGMENTS 

We thank Albert Bosnia for useful discussions on the data. 
This work was partially supported by grant ANR-06-BLAN- 
0172, by the Russian Foundation for Basic Research (grants 
08-02-00361-a and 09-02-00968-a) and by a grant from the 
President of the Russian Federation for support of Leading 
Scientific Schools (grant NSh-1318.2008.2). 



REFERENCES 

Athanassoula E., Vozikis Ch., 1999, in Barnes J. E., 
Sanders D.B., eds, Proc. IAU Symp. 186, Galaxy inter- 
actions at low and high redshift, Kluwer Academic Pub- 
lishers, Dordrecht, p. 145 

Athanassoula E., 2005a, in Szczerba R., Stasihska G. & 
Gorny S. K., eds, AIP Conference Proceedings 804, Plan- 
etary Nebulae as Astronomical Tools, Melville, New York, 
p. 333 

Athanassoula E., 2005b, MNRAS, 358, 1477 

Binney J., Tremaine S., 2008, Galactic Dynamics: Second 
Edition. Princeton Univ. Press, Princeton 

Bosma A., 2004, in Ryder S.D., Pisano D.J., Walker M.A., 



Freeman K.C., Proc. IAU Symp. 220, Dark matter in 
Galaxies, Astron. Soc. Pac, San Francisco, p. 39 

Coccato L., Gerhard O., Arnaboldi M., Das P., Douglas N. 
G., Kuijken K., Merrifield M. R., Napolitano N. R., Noor- 
dermeer E., Romanowsky A.J., Capaccioli M., Cortesi A., 
de Lorenzi F., Freeman K. C, 2009, MNRAS, 394, 1249 

de Lorenzi F., Debattista V. P., Gerhard O., Sambhus N., 
2007, MNRAS, 376, 71 

de Lorenzi F.; Gerhard O., Saglia R. P., Sambhus N., De- 
battista V. P., Pannella M., Mendez R. H., 2008, MNRAS, 
385, 1729 

de Lorenzi F., Gerhard O., Coccato L., Arnaboldi M., Ca- 
paccioli M., Douglas N. G., Freeman K. C, Kuijken K., 
Merrifield M. R., Napolitano N. R., Noordermeer E., Ro- 
manowsky A. J., Debattista V. P., 2009, MNRAS, 395, 
76 

Dehnen W., 2000, ApJ, 536, 39 

Dehnen W., 2000, AJ, 119, 800 

Dehnen W., 2002, Journal of Computational Physics, 179, 
27 

Dekel A., Stoehr F., Mamon G. A., Cox T. J., Novak G. 
S., Primack J. R., 2005, Nature, 437, 707 

Douglas N. G., Napolitano N. R., Romanowsky A. J., Coc- 
cato L., Kuijken K., Merrifield M. R., Arnaboldi M., Ger- 
hard O., Freeman K. C, Merritt H. R., Noordermeer E., 
Capaccioli M., 2007, ApJ, 664, 257 

Gerhard O.E., 1993, MNRAS, 265, 213 

Goudfrooij P., Hansen L., Jorgensen H. E., Norgaard- 
Nielsen H. U., de Jong T., van den Hoek L. B., 1994, 
A&AS, 104, 179 

Joanes D.N., Gill C.A., 1998, The Statistician, 47, 183 

Lauer T. R., Faber S. M., Gebhardt K. et al., 2005, AJ, 
129, 2138 

Machado R. E. G. and Athanassoula E., 2010, MNRAS, 
submitted 

Mendez R. H., Riffeser A., Kudritzki R.-P., Matthias M., 
Freeman K. C, Arnaboldi M., Capaccioli M., and Gerhard 
O. E., 2001, ApJ, 563, 135 

Napolitano N. R., Romanowsky A. J., Coccato L., Capac- 
cioli M., Douglas N. G., Noordermeer E., Gerhard O., 
Arnaboldi M., de Lorenzi F., Kuijken K., Merrifield M. 
R., O'Sullivan E., Cortesi A., Das P., Freeman, K. C, 
2009, MNRAS, 393, 329 (N09) 

Navarro J.F., Frenk C.S., White S.D.M., 1996, ApJ, 462, 
563 

Navarro J.F., Frenk C.S., White S.D.M., 1997, ApJ, 490, 
493 

Kenney J. F., Keeping, E. S. "Mathematics of Statistics", 
Pt. 2, 2nd ed. Princeton, NJ: Van Nostrand, 1951. 

Rodionov S.A., Orlov V.V., 2008, MNRAS, 385, 200 

Rodionov S. A., Athanassoula, E., Sotnikova, N. Ya. 2009, 
MNRAS, 392, 904 (RAS09) 

Romanowsky A. J., Douglas N. G., Arnaboldi M., Kuijken 
K., Merrifield M. R., Napolitano N. R., Capaccioli M., 
Freeman K. C, 2003, Science, 301, 1696 

Sersic J. L., 1968, Atlas de galaxias australes. Observatorio 
Astronomico, Cordoba, Argentina 

Syer D., Tremaine S., 1996, MNRAS, 282, 223 

van der Marel R. P., Franx M., 1993, ApJ, 407, 525 

Weil M. L., Hernquist L., 1994, ApJ, 431, 79 

Weil M. L., and Hernquist L., 1996, ApJ, 460, 101 



Dynamical models of the elliptical galaxy NGC 44 $4- 17 



APPENDIX A: CALCULATION OF 
GAUSS-HERMITE MOMENTS IN THE CASE 
OF A-BODY SYSTEMS 

Let us now describe how we calculate the Gauss- 
Hermite moments in the cas e of A-body systems . Both 
Ivan der Marel fc Franxl l| 19931 ) and IGerhardl (|l993l ). inde- 
pendently, discuss the use of Gauss-Hermite moments to 
measure the deviations of the observed velocity distribu- 
tion from a Gaussian, but with different normalization of 
the Hermite po lynomials. Here we use the normalization of 
IGerhardl j 19931 ). 

The Hermite polynomials are defined as 

H n (x) = (-ire* 2 -^(e-* 2 ) (Al) 

The sequence of Hermite polynomials satisfies the recursion 



H n +i — 2xH n 



2nH n 



and the first three Hermite polynomials are 



Ho 



1, Hi = 2x, H 2 = 4aT 



The set of functions defined as 

u n (x) = (2 n+1 n\7vy 1/2 H n (x)e X p(-x 2 /2) . 

obey the orthogonality relation 

+ oo en 

U n (x)u m (x)dx = 



27T 1 / 2 ' 



(A2) 



(A3) 



(A4) 



(A5) 



Thus this set of functions is a complete orthogonal system. 
The Gauss-Hermite moments for some function l(v) are 
defined as 

h„ = 27r 1/2 7^ 1 / l(v)u n (w)dv , w = (v- V h )/a h (A6) 

J — oo 

where Vh, (?h 7^ and 7/1 7^ are free parameters. The 
function l(v) can be approximately calculated by means of 
a truncated Gauss-Hermite series 



l(v) « — - S~]hjUj(i 



<7h 



(A7) 



whe re Nh is the number of t erms used (see IGerhardl (|l993T ) 
and Ivan der Marel fc Franxl (1993) for more details). 

In the case of an A-body system we need to solve the 
following problem. We have the set of values Vi, V2 ■■ ■ Vn, 
and we need to calculate the Gauss-Hermite moments of the 
distribution function l(v) of these values. For example, if we 
need to calculate Gauss-Hermite moments for the line-of- 
sight velocity distribution in some area of an iV-body sys- 
tem then Vi is the line-of-sight velocity of each particle in 
this area. In this case, the Gauss-Hermite moments for the 
function l(v) can approximately be calculated as 



2 1/2 n 



i) , Wi = (vi - V h )/a h . 



(A8) 



We use this equation as the definition of Gauss-Hermite mo- 
ments for our discrete case. 

To compare Gauss-Hermite moments with the observa- 
tions, we need to choose the free parameters Vh, o~h and jh as 
observer would do. These free para meters should be chosen 
so as to give ho = 1, hi — hi = (|van der Marel fc Franxl 



Il993h . From |[A3]). (1A4)) , (|A8)l for n = 1 and the condition 
hi — 0, we have 

E^Ii v i exp(-iu I 2 /2) 



hi = <£> V h = 



^N 



Er =1 exp(-» 2 /2) 



(A9) 



From (|A3|) . (|A4|) (IA8|) for n — 2 and the condition h 2 = 0, 
we have 

h 2 - <^ o h - 2 — w — . (A10) 

Ei=i exp(-w 2 /2) 

We note that Wi = (vt — Vh)/<7h, so ()A9|I and (|A10[) cannot 
be solved directly. We solve (|A9|l and (|A10[) together, by 
means of iterations. Initially we set Vh equal to the mean 
value of Vi, (Jh equal to the standard deviation of Vi and 
7ft = 1. A single iteration is as follows 



v; 



{new) 



E^Li^ ex p(-<7 2 ) 

Ef =1 exp(- W 2 /2) 



(new) 



= 2 



,- ^-^) 2 e X p(~^/2 ) V:J !XMi 
E l J I 1 exp(-» 2 /2) 



After finding the appropriate Vh and an values, we can easily 
calculate the last free parameter as fh ew — ho- 



APPENDIX B: CALCULATION OF SKEWNESS 
AND KURTOSIS FOR A SMALL SAMPLE 

In the outer parts of the galaxy we have information about 
the line-of-sight velocity distribution only from PNe obser- 
vations. Our task is to compare the velocity distribution of 
the PNe in some area on the sky plane with the velocity dis- 
tribution in a constructed A r -body model (see Fig. [T} . More 
precisely our task can be formulated as follows. We have 
two random samples which we denote here as "A" and "B" . 
Sample "A" with size n a consist of the line-of-sight veloci- 
ties of the PNe in the selected area. Sample "B" with size n;, 
consist of the line-of-sight velocities of model particles in the 
same area. We need to assess the probability that these two 
samples were generated from the same distribution function. 
This can be achieved by comparing moments calculated for 
these two samples. 

Here we will discuss the comparison of the skewness 
and kurtosis calculated for these two samples. In our case 
the number of planetary nebulae in a given area (the size of 
the sample A) is small, not exceeding 18 (Table [TJ. On the 
other hand the size of sample B, i.e. the number of parti- 
cles in given area of constructed model is rather large, and 
is about 10 4 . So we need to compare high-order moments 
calculated for small and for large samples. Of course, for 
a sample as small as A the uncertainties of the estimators 
of the higher-order moments are rather large. Moreover, we 
cannot calculate these uncertainties without making an as- 
sumption about the distribution function. If we assume, for 
example, a Gau ssian distribution t hen we can calculate these 
uncertainties (|Joanes fc Gilllll998l ). but if we do not assume 
a priori any distribution function then the uncertainties are 
formally infinite. Furthermore, in the general case all com- 
monly use d estimators of the s ample skewness and kurtosis 
are biased (|Joanes fc Gillll 19981 ) . For small samples this bias 
can be very large, esp e cially for the kurtosis (see tables 2 
and 3 in ljoanes fc Gilll (J1998)). Consequently, for sample A 
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the bias can be rather large. The bias depends on the size of 
the sample and for sample B (which is rather large) the bias 
is negligibly small. If both of our samples are generated from 
the same distribution function then the expected values of 
kurtosis (and skewness) for these two samples can differ sig- 
nificantly. We can solve this problem reducing sample "B" 
to a smaller sample size. 

Let us now describe how we calculate kurtosis. The 
skewness is calculated in the same way. Let us have some 
estimator of sample kurtosis. For sample A we simply cal- 
culate value of kurtosis K a using the chosen estimator. For 
sample B we calculate the value of the kurtosis "reducing" 
the sample size to that of sample A, i.e. n a . We randomly get 
n a members from sample B and calculate for this sub-sample 
the kurtosis k\. We repeat this Nk times and calculate fci, 
k 2 ... k Nk (in this article N k = 100 000 ). We denote the 
mean value of ki as K™ a and the standard deviation of ki 
as (T™ a . Let fb{x) be the distribution function corresponding 
to sample B. For the chosen kurtosis estimator we can con- 
struct the distribution function gb(x) of the sample kurtosis 
for samples with size n a generated from fb(x). We note that 
K™ a is approximately equal to the expected value of gb(x) 
and (7^ is approximately equal to the standard deviation of 
gb(x). For example, if the sample A was also generated from 
the distribution function fb(x), then K^ a would be the ex- 
pected value and a^ a the standard deviation of the kurtosis 
of sample A. 

We use K a as measure of the kurtosis of sample A, K % 
as measure of the kurtosis of sample B, and a^ as measure 
of the standard deviation. We note that t h e mai n value for 
us is r = \K a - AT|/q"°. |joanes fc Gilll (|l998T l discussed 
three different estimators of the kurtosis. For our analysis, 
however, it doesn't matter which of the three we use be- 
cause it can be proven that the value of r is the same for 
all of them. To calculate the sample skewness and the sam- 
ple kurtosis w e use t he estimators Gi and G2, discussed in 
Ijoanes fc Gill (1998). These estimators are unbiased in the 
case of a normal distribution. 

As we noted above, we cannot calculate uncertainties 
of estimators of kurtosis or skewness without making an as- 
sumption about the distribution function. The value a c b l is 
the standard deviation of the estimation of the kurtosis or 
the skewness for sample A (observation) assuming that it s 
from the same distribution function as sample B (model). So 
this value depends on the model. We also can calculate the 
standard deviations <rJJ of the estimation of kurtosis or skew- 
ness for sa mple A assuming th at its distribution function is 
Gaussian IjJoanes fc Gilll Il998l ). For our models in general 
value of (7^ is less that of a^. The difference, however, is 
not so strong. The maximum difference is cr^/cr^ w 1.3 and 
a b l G n ~ 1-4 for the skewness and the kurtosis, respectively. 



